PROGRAM convert(input,output);

{convert decagonal diffraction data to approximant data}
{splitted from mrq.p}
{Version 7: Chemnitz Oct 96}

{
 WARNING 1.0 : for 2D RECTANGULAR (3D orthorhombic) unit cells only!!!
 WARNING 1.1 : be careful to provide 5d coordinates in proper format
 X-axis: [2*Fx(0), Fx(-1), -Fx(1), -Fx(1), Fx(-1)]
 Y-axis: [      0,  Fy(1),  Fy(0), -Fy(0), -Fy(1)]
 
 INPUT: TILING (file with 5d coordinates of tiling objects)
 INPUT: expe data
 INPUT: mat (transform matrix for 5d expe vectors -
              comments therein)
 OUTPUT: dat.hkl
 
}

CONST
   appnum = 5; {number of the independent approx. HK's making deca orbit}
   edata  = 4000; {number of datapoints}
   ndatap = edata;
   tsym	  = 40;

{LOCAL STUFF}
TYPE
   MYREAL   = single;
   mystring = string[100];
   v5	    = array[1..5] of integer;
   r5	    = array[1..5] of MYREAL;
   mat5x5   = array[1..5] of v5;
   matr5x5  = array[1..5] of r5;
   v2	    = array[1..2] of MYREAL;
   v3	    = array[1..3] of MYREAL;
   i2	    = array[1..2] of integer;
   i3	    = array[1..3] of integer;   
   i10	    = array[1..10] of integer;
{GAUSSJ, LFIT, MRQMIN, MRQCOF, COVRST...}
   
VAR

   txt : text;
   tiling_name,diff_data,ch : mystring;
   
   ndat,o,h,k,l,i,j,maxtab,h0,h1,k0,k1:integer;
   nsr10,ns2,multi:integer;
   
   hx,kx,lx,emax,imax,xx,eps,zero,qq,pi,por,tau,eta,sig,sh,
   nsteu,arg,s_plane,s_per:myreal;
   
   b:V5;
   v,phason:V2;
   iv:I2;
   bb:ARRAY[1..2] of v5;
   bs:ARRAY[1..2] of r5;
   ba,br,bap:ARRAY[1..2] of v2;
   e,es,ep,esp,ed : ARRAY[1..5] of v2;
   edip,edi : ARRAY[1..5] of i2;
   sr10 : ARRAY[1..tsym] of v5;
   s2 : ARRAY[1..tsym] of i2;
   hke : ARRAY[1..edata,1..tsym] of i2;
   le,napp : ARRAY[1..edata] of INTEGER;
   be : ARRAY[1..edata] of v5;
   sm:matr5x5;
   what,recmirror: integer;
   c,s:r5;
   muco:real;

procedure tra2(x : v2;Var y:v2);
begin
   y[1]:=(x[1]*ba[2,2]-x[2]*ba[2,1])/(ba[1,1]*ba[2,2]-ba[1,2]*ba[2,1]);
   y[2]:=(x[2]*ba[1,1]-x[1]*ba[1,2])/(ba[1,1]*ba[2,2]-ba[1,2]*ba[2,1]);
end; { tra }
procedure tra2r(x : v2;Var y:v2);
begin
   y[1]:=(x[1]*br[2,2]-x[2]*br[2,1])/(br[1,1]*br[2,2]-br[1,2]*br[2,1]);
   y[2]:=(x[2]*br[1,1]-x[1]*br[1,2])/(br[1,1]*br[2,2]-br[1,2]*br[2,1]);
end; { tra }

function absv(x	: v2):MYREAL;
begin
   absv:=sqrt(x[1]*x[1]+x[2]*x[2]);
end;
   
procedure union2(x : i2);
var o,i,j	    : INTEGER;
begin
   o:=0;
   repeat
      o:=o+1;
      i:=0; for j:=1 to 2 do if (abs(s2[o,j])=abs(x[j])) then i:=i+1;
   until (o=ns2) or (i=2);
   if (o=ns2) and (i<>2) then
   begin
      ns2:=ns2+1;
      s2[ns2]:=x;
   end;
end;
procedure union2all(x : i2);
var o,i,j	    : INTEGER;
begin
   o:=0;
   repeat
      o:=o+1;
      i:=0; for j:=1 to 2 do if (s2[o,j]=x[j]) then i:=i+1;
   until (o=ns2) or (i=2);
   if (o=ns2) and (i<>2) then
   begin
      ns2:=ns2+1;
      s2[ns2]:=x;
   end;
end;
procedure union5(x : v5);
var o,i,j	    : INTEGER;
begin
   o:=0;
   repeat
      o:=o+1;
      i:=0; for j:=1 to 5 do if (sr10[o,j]=x[j]) then i:=i+1;
   until (o=nsr10) or (i=5);
   if (o=nsr10) and (i<>5) then
   begin
      nsr10:=nsr10+1;
      if (nsr10>20) then writeln('hey! nsr10 = ',nsr10:1);
      sr10[nsr10]:=x;
   end;
end;
procedure mir_per(x : v5);
var
   a,b		    : v5;
   o		    :  INTEGER;
begin
   a:=x;
   o:=0;
   repeat
      o:=o+1;
      union5(a);
      b:=a;
      a[1]:=-b[1]; a[2]:=-b[5]; a[3]:=-b[4]; a[4]:=-b[3]; a[5]:=-b[2];
   until o=2;
end; { mir_per }
procedure mir_par(x : v5);
var
   a,b : v5;
   o   : INTEGER;
begin
   a:=x;
   o:=0;
   repeat
      o:=o+1;
      mir_per(a);
      b:=a;
      a[1]:=b[1]; a[2]:=b[5]; a[3]:=b[4]; a[4]:=b[3]; a[5]:=b[2];
   until o=2;
end;
procedure decam5(x :  v5);
var
   a,b : v5;
   o   : INTEGER;
begin
   a:=x;
   {5-fold rotation}
   o:=0;
   repeat
      o:=o+1;
      if (recmirror>0) then mir_par(a) else union5(a);
      b:=a;
      a[1]:=b[5]; a[2]:=b[1]; a[3]:=b[2]; a[4]:=b[3]; a[5]:=b[4];
   until o=5;
end; { decam }
procedure decam10(x : v5);
var
   a,b : v5;
   o   : INTEGER;
begin
   a:=x;
   {5-fold rotation}
   o:=0;
   repeat
      o:=o+1;
      if (recmirror>0) then mir_par(a) else union5(a);
      b:=a;
      a[1]:=-b[4]; a[2]:=-b[5]; a[3]:=-b[1]; a[4]:=-b[2]; a[5]:=-b[3];
   until o=10;
end; { decam }
procedure read_expe;
var o,i,j,b5		   : INTEGER;
   y,pi2,qrms		   : MYREAL;
   ww,vv,dv,vp		   : V2;
   txt			   : text;
   siginp,nde		   : integer ;
   ie,er,ic,erc,q,dq,ex,qp : ARRAY [1..ndatap] OF MYREAL;
begin
   pi2:=2.0*Pi;
   reset(txt,diff_data);
   readln(txt,siginp);
   readln(txt); readln(txt);
   emax:=0.0;
   nde:=0;
   qrms:=0;
   repeat
      ndat:=ndat+1;
      for o:=1 to 4 do read(txt,b[o]); readln(txt,b5,ie[ndat],y);
      if (siginp=1) then er[ndat]:=1/sqrt(y) else er[ndat]:=y;
      b[5]:=0;
      be[ndat]:=b;
      be[ndat,5]:=b5; {replace 5th deca coord by periodic}
      if (emax<ie[ndat]) then emax:=ie[ndat];

      {approximant peak position}
      for o:=1 to 2 do iv[o]:=0;
      for o:=1 to 2 do
         for j:=1 to 5 do iv[o]:=iv[o]+edi[j,o]*b[j];
      hke[ndat,1,1]:=iv[1];
      hke[ndat,1,2]:=iv[2];
      {hke[ndat,1,3]:=b5;}
      {vv[1]:=hx*iv[1]/s_plane/pi2;
      vv[2]:=kx*iv[2]/s_plane/pi2;}
      vv[1]:=(iv[1]*br[1,1]+iv[2]*br[2,1])/s_plane;
      vv[2]:=(iv[1]*br[1,2]+iv[2]*br[2,2])/s_plane;
      xx:=b5*1.0/s_per;

      {quasicrystal peak position, par and perp}
      for o:=1 to 2 do begin v[o]:=0.0; vp[o]:=0.0; end;
      for o:=1 to 2 do
      begin
	 for j:=1 to 5 do v[o]:=v[o]+es[j,o]*b[j]/s_plane;
	 for j:=1 to 5 do vp[o]:=vp[o]+esp[j,o]*b[j]/s_plane;
      end;
      le[ndat]:=b5*multi;
      q[ndat]:=sqrt(v[1]*v[1]+v[2]*v[2]+xx*xx);

      {for o:=1 to 5 do write(b[o]:5);
      writeln(2.0*pi*v[1]:8:3,2.0*pi*v[2]:8:3,2.0*pi*q[ndat]:8:3);}
{      decam10(b); writeln(nsr10:1); for i:=1 to nsr10 do begin
	 for o:=1 to 5 do write(sr10[i,o]:5); writeln;
      end;}

      qp[ndat]:=absv(vp);

      {displacement}
      for o:=1 to 2 do dv[o]:=vv[o]-v[o];
      dq[ndat]:=absv(dv);
      qrms:=qrms+dq[ndat]*dq[ndat];

      {generating all independent approximant peaks}
      nsr10:=0;
      sr10[1]:=b; nsr10:=1;
      decam10(b);
      napp[ndat]:=nsr10;
      {writeln(ndat:1,' ',nsr10:1);}
      for i:=1 to nsr10 do
      begin
	 for o:=1 to 2 do iv[o]:=0;
	 for o:=1 to 2 do
	    for j:=1 to 5 do iv[o]:=iv[o]+edi[j,o]*sr10[i,j];
	 hke[ndat,i]:=iv;
	 {
	 if(b5=0) then begin
	 write('TEST ',i:1,'  ');
	 for o:=1 to 5 do write(sr10[i,o]:4);
	 write('  ');
	 for o:=1 to 2 do write(iv[o]:4);
	 ww[1]:=(iv[1]*br[1,1]+iv[2]*br[2,1])/s_plane;
	 ww[2]:=(iv[1]*br[1,2]+iv[2]*br[2,2])/s_plane;
	 writeln('  ',ww[1]:8:4,ww[2]:8:4,'  |q|=',sqrt(ww[1]*ww[1]+ww[2]*ww[2]):1:5,' ndat=',ndat:1,' ii=',ie[ndat]:1:5);
	 end;
	 }
      end;
      napp[ndat]:=nsr10;
      if (what=0) then napp[ndat]:=1;
      if (what>0) then begin
	 s2[1]:=hke[ndat,1]; ns2:=1;
	 if(what=1) then for i:=1 to napp[ndat] do union2(hke[ndat,i]);
	 if(what=2) then for i:=1 to napp[ndat] do union2all(hke[ndat,i]);
	 napp[ndat]:=ns2;
	 for i:=1 to ns2 do hke[ndat,i]:=s2[i];
      end;
   until eof(txt);
   close(txt);

   {get RMS Q_displacement approximant vs quasicrystal}
   qrms:=sqrt(qrms/ndat);
   
   {writeln('Experimental data: ',ndat:1,' reflections.');}
   {normalize to the brightest peak}
   for i:=1 to ndat do
   begin ie[i]:=ie[i]/emax; er[i]:=er[i]/emax; end;
   writeln(ndat:1,'  ',br[1,1]/s_plane:1:7,' ',br[1,2]/s_plane:1:7,' ',br[2,1]/s_plane:1:7,' ',br[2,2]/s_plane:1:7,' ',1.0/s_per/multi:1:7,'  # |F| wt=1/sig^2 hkl1 h5 |Q| |Q_perp| |Q-Q_app| napp hkls');
   for i:=1 to ndat do begin
      write(ie[i]:1:6,'  ',er[i]:1:6,'  ');
      write(hke[i,1,1]:1,' ',hke[i,1,2]:1,' ',le[i]:1,'  ');
      for o:=1 to 5 do write(be[i,o]:1,' ');
      write(q[i]:7:4,qp[i]:7:4,dq[i]:7:4,' ',(napp[i]-1):1,'  ');
      for j:=2 to napp[i] do begin
	 write(hke[i,j,1]:1,' ',hke[i,j,2]:1,' ',le[i]:1,' ');
	 write(' ');
      end;
      writeln;
   end;
end;
procedure READ_TILING_BASIS;
var o,i,j,k,l: integer;
   fnplus,fnminus,fn0,fmplus,fmminus,fm0: myreal;
begin
  {read input file for the tiling}
   reset(txt,tiling_name);   
   for i:=1 to 2 do
   begin
      for o:=1 to 2 do begin ba[i,o]:=0.0; bap[i,o]:=0.0; end;
{5D PROJECTION: assuming first basis vector is along X and second along Y}
      for o:=1 to 5 do read(txt,b[o]); readln(txt);
      for j:=1 to 2 do
	 for o:=1 to 5 do
	 begin
	    ba[i,j]:=ba[i,j]+b[o]*e[o,j];
	    bap[i,j]:=bap[i,j]+b[o]*ep[o,j];
	 end;
      k:=b[1]; for o:=1 to 5 do b[o]:=b[o]-k;
      bb[i]:=b;
      phason[i]:=absv(bap[i]);
   end;
   readln(txt,o,multi);
   {reciprocal basis vectors}
   eps:=0.000001;
   por:=abs(ba[1,1]*ba[2,2]-ba[2,1]*ba[1,2]);
   br[1,1]:=ba[2,2]/por;
   br[2,2]:=ba[1,1]/por;
   br[1,2]:=-ba[2,1]/por;
   br[2,1]:=-ba[1,2]/por;
   por:=br[1,1]*br[2,1]+br[1,2]*br[2,2];
   {
   if abs(por)>eps then
   begin
      writeln('Non-rectangular unit cell. Sorry, halting.');
      halt;
   end;
   if (((abs(br[1,1])>eps) and (abs(br[1,2])>eps)) or
      ((abs(br[2,1])>eps) and (abs(br[2,2])>eps))) then
   begin
      writeln('Improperly oriented unit cell!');
      writeln('Sorry, halting.');
      halt;
   end;
   }
   if (abs(br[1,1])>eps) then hx:=br[1,1];
   if (abs(br[1,2])>eps) then kx:=br[1,2];
   if (abs(br[2,1])>eps) then hx:=br[2,1];
   if (abs(br[2,2])>eps) then kx:=br[2,2];

   for i:=1 to 5 do
   begin
      {edi[i,1]:=round(es[i,1]/hx); edi[i,2]:=round(es[i,2]/kx);
      writeln('edi  ',i:1,edi[i,1]:5,edi[i,2]:5);}
      tra2r(es[i],v);
      edi[i,1]:=round(v[1]); edi[i,2]:=round(v[2]);
      {writeln('edi2 ',i:1,edi[i,1]:5,edi[i,2]:5);}
      {writeln('es ',i:1,es[i,1]:9:5,es[i,2]:9:5);}
   end;

   hx:=2.0*pi*hx;
   kx:=2.0*pi*kx;
   lx:=2.0*pi;
   close(txt);
end; { READ_TILING }

BEGIN
   {what==0: export only 1 hkl datapoint for each expe datapoint}
   {what==1: export symmetry-inequivalent hkl's}
   {what==2: export all hkl's}
   readln(what);
   readln(recmirror);
   readln(diff_data);  {name of the datafile}
   readln(tiling_name); {name of the approximant file - read basis}
   
   pi:=4.0*arctan(1.0);
   tau:=(1+sqrt(5.0))/2.0;
   eta:=sqrt(tau+2.0);
   nsteu:=2.0*tau/5.0*tau*tau; {Unused. Just to remember...}
      
   {pentagonal basis (normalized to the edge length of the DPH tiling)}
   for i:=1 to 5 do
   begin
      e[i,1]:=cos(2.0*pi*(i-1.0)/5.0); 
      e[i,2]:=sin(2.0*pi*(i-1.0)/5.0);
   end;

   ep[1]:=e[1];
   ep[4]:=e[2];
   ep[2]:=e[3];
   ep[5]:=e[4];
   ep[3]:=e[5];


   {get unit cell}
   reset(txt,'INPUT/PPOT/uc.inp');
   readln(txt,s_plane,s_per);
   close(txt);
   {writeln('Lattice parameters: ',s_plane:1:3,' ',s_per:1:3);}
   
   {get matrix transforming lattice translations into DIF_DATA convention}
   reset(txt,'INPUT/DFIT/mat.hkl');
   for i:=1 to 5 do
   begin
      for j:=1 to 5 do read(txt,sm[i,j]); readln(txt);
      {for j:=1 to 5 do sm[i,j]:=sm[i,j]*2/5/tau;}
   end;
   close(txt);

   {Diffraction data convention (see comments in *mat}
   for j:=1 to 5 do
      for o:=1 to 2 do es[j,o]:=0.0;
   for j:=1 to 5 do {cycle through 5 star vectors}
      for i:=1 to 5 do
      begin
	 es[j,1]:=es[j,1]+sm[j,i]*e[i,1];
	 es[j,2]:=es[j,2]+sm[j,i]*e[i,2];
	 {por:=es[j,1]; es[j,1]:=es[j,2]; es[j,2]:=por;}
      end;
   for j:=1 to 5 do
   begin
      por:=absv(es[j]);
      {writeln('es[',j:1,'] : ',es[j,1]:6:1,es[j,2]:6:1,' length ',por:1:5);}
   end;
   {steurer}
   
   for i:=1 to 5 do c[i]:=cos(2*pi*(i-1)/5);
   for i:=1 to 5 do s[i]:=sin(2*pi*(i-1)/5);
   
   {muco:=2/5*tau*tau; 
    for i:=1 to 5 do begin
    for j:=1 to 2 do begin
    es[i,j]:=e[i,j]*muco;
    esp[i,j]:=esp[i,j];
      end;
   end;}
   
   esp[1]:=es[1];
   esp[4]:=es[2];
   esp[2]:=es[3];
   esp[5]:=es[4];
   esp[3]:=es[5];
    
   read_tiling_basis;

   {EXPERIMENTAL DATA}
   read_expe;

END.


